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ABSTRACT 

In order to reproduce the statistical properties of the observed exoplanets, population synthesis models have 
shown that the migration of protoplanets should be significantly slowed down, and that processes stalling mi- 
gration should be at work. Much current theoretical efforts have thus been dedicated to find physical effects that 
slow down, halt or even reverse migration. Many of these studies rely on the horseshoe drag, whose long term- 
evolution (saturated or not) is intimately related to the disk viscosity in laminar disk models. We investigate 
how the horseshoe drag exerted on a low-mass planet is altered by a more realistic treatment of the turbulence 
in protoplanetary disks. Two-dimensional hydrodynamic simulations are performed with a turbulence model 
that reproduces the main turbulence properties of three-dimensional magnetohydrodynamic calculations. We 
find that the horseshoe drag can remain unsaturated on the long term, depending on the turbulence strength. We 
show that the desaturation of the horseshoe drag by turbulence can be modeled by vortensity diffusion across 
the time-averaged planet's horseshoe region. At low-turbulence, the running-time averaged torque is in good 
agreement with the total torque obtained for an equivalent laminar model, with a similar vortensity's diffusion 
coefficient. At high-turbulence, differences arise due to the time-evolution of the averaged density profile with 
turbulence. 

Subject headings: accretion, accretion disks — hydrodynamics — turbulence — methods: numerical — plan- 
etary systems: formation — planetary systems: protoplanetary disks 



L INTRODUCTION 

The discovery of exoplanets has provided an exciting op- 
portunity to test our theories of planet formation and evolu- 
tion. T he detection of the first hot Jupiter by Mayor & Ouelo3 
(Il995h made planetary migration become a central ingredient 
of these theories (Lin eLal. 1996). Planetary migration occurs 
from the early stages of planet formation, as planets tidally 
interact with the protoplanetary disk where they are embed- 
ded in (iGoldreich & Tremaindll980h . Planets experience a 
torque from the disk that alters their eccentricity and semi- 
major axis. The variation of the planets semi-major axis under 
the disk torque is referred to as planetary migration. 

Depending on the disk and planet properties, three types 
of migration are usually distinguished. Type I migration 
applies to low-mass planets, typically up to a few Earth 
masses if the central object has a solar mass. In isother- 
mal disks, type I migration is directed inwards and is much 
faster than both the f ormation of giant planets cores and the 
disk evaporation (e.g. lWardllI997l:lTanaka et al.ll2002l) . To re- 
produce the mass-distance distribution of the observed exo- 
planets, population synthesis models have shown that type 
I migration should be significantly slowed down (llda & LinI 
|2008a; Mordasini et al. 2009), and that processes stalling mi- 
gration should be at work, in particular ne ar the snow line 
dlda & LinI l2008bt ISchlaufman eTall l2009l) . Some mech- 
anisms recently proposed to slow down, halt or even re- 
verse migration are reviewed below. Type II migration con- 
cerns planets that are massive enough to open up a clean 
gap around their orbit (typically Jupiter-mass planets). In 
this regime, planets dri ft inwards at the disk viscous rate 
jLin & Papaloizoull986l) . Intermediate-mass planets building 
up a partia l gap can u ndergo a runaway migration in massive 
disks (Mas set & Papa loizou 2003). This is usually referred to 

Electronic address: [clement.baruteau@ucolick.org; lin@ ucolick.org | 



as type III migration. Throughout this paper, we will focus on 
type I migration. 

1. 1. Type I migration 

The torque exerted on a low-mass planet comprises the 
differential Lindblad torque and the horseshoe drag (e.g. 
iPaardekooper & Papaloizoiill2009al) . The former results from 
the exchange of angular momentum between the planet and 
the circulating^ fluid elements located at Lindblad resonances 
( IGoldreich & Tremaind [l979h^ The differential Lindblad 



torque is generally negative (IWardI 19971). and is stationary 
after a few dynamical timescales (iMeyer-Vernet & Sicardyl 
119871) . Its value is however particularly sensitive to any mech- 
anism altering th e location of the resonances, such as the 
disk self-gravity dPierens & Hurel 120051: iBaruteau & Massej 
2008b). 

The horseshoe drag accounts for the exchange of angu- 
lar momentum between the planet and the li bratin g ' fluid 
elements in the planet 's horseshoe region (IWardr il991). 
iCasoh & Massej (120091) have clarified that the horseshoe 
drag actually corresponds to the corotation torque, namely 
the torque exerted by the corotation region^ on the planet. 
The horseshoe drag is a non-li near phenomenon th at can- 
not be captured by linear theory ( IPaardekooper & Papaloizoul 
l2009ah . except at early times, before fluid streamUnes get 
fully adjusted to the planet's introduction in a timescale com- 
parable to the horseshoe U-turn time (IBaruteau & Massed 
2008a). Contrary to the differential Lindblad torque, the 
horseshoe drag can be either positive or negative, and it is 

' In the frame corotating with the planet. 

^ As pointed out by Casoli & Masse| <2009l) . the corotation and horseshoe 
regions generally differ. While a linear analysis indicates that the corotation 
region has a radial width of approximately a pressure scaleheight H about the 
corotation radius, the width of the horseshoe region vanishes as the planet's 
mass tends to zero. 
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generally not a stationary quantity. In isothermal disks, its 
time dependence is related to the time evolution of the disk 
vortensity (the vorticity to surface density ratio) in the horse- 
shoe region. The maximum value of the horseshoe drag 
scales with the (opposite of the) unperturbed vortensity gra- 
dient across the horseshoe region (Ward 1992; Masset 2001'). 
This maximum value is called the fully unsaturated horseshoe 
drag. It is positive for surface density profiles locally shal- 
lower than r^^l^. In the particular case of a planet fixed in 
an inviscid disk, the horseshoe region is closed and vorten- 
sity is conserved along horseshoe streamlines. As vorten- 
sity is progressively stirred up in the horseshoe region, the 
horseshoe drag oscillates with time with a decreasing ampli- 
tude. The horseshoe drag ultimately cancels out as vorten- 
sity gets uniformly distributed after several libration times 
(|Balmforth & Korvcanskv 2001). This is known as the horse- 
shoe drag saturation. 

For radiative disks in the adiabatic limit, the horseshoe 
drag has an additional contribution, whose maximum value 
scales with the (opposite of the) unperturbed entropy gradi- 
ent across the horseshoe regio n (Baruteau & Masset 2008a; 
iPaardekooper & Mellemal2008i;lMasset & Casoli 2009). This 
additional contribution is called entropy-related horseshoe 
drag. The contribution already existing in isothermal disks, 
and which features the vortensity gradient, is similarly de- 
noted as vortensity-related horseshoe drag for future refer- 
ence. The entropy-related horseshoe drag is positive for neg- 
ative entropy gradients. It can therefore slo w down, halt or 
even reverse (IPaardekooper & Mellemal 12006 ) type I migra- 
tion compared to isothermal disk models, as long as saturation 
does not occur. 

1.2. Horseshoe drag desatumtion 

As described above, the horseshoe drag in isothermal disks 
saturates as vortensity is strictly advected along horseshoe 
streamlines. Viscosity, which acts as a source t erm in the 
vortensity equation (e.g. lOgilvie &Lubowll2006h . can sus- 
tain a non-zero vortensity gradient across the horseshoe re- 
gion. If so, the vortensity-relat ed horsesh oe drag would attain 
a steady-state value (Masset 2001', '2002), which arises from 
a net exchange of angular momentum between the horseshoe 
region and the rest of the disk. In particular, the torque re- 
mains unsaturated when the viscous diffusion time across the 
horseshoe r egion is somewhat lar ger than the horseshoe U- 
turn time (Baruteau & Masset 2008a), but smaller than the li- 
bration time (e.g. Ward 1991; Masset 2001). At large viscosi- 
ties, the horseshoe drag decreases with increasing viscosity, 
as viscous diffusion tends to impose the initial vortensity pro- 
file. This behavior is known as the cut-off of the horseshoe 
drag (Masset"2001|). For strong enough viscosities, the horse- 
shoe drag coin cides with the line ar estimate of the corotation 
torque (IPaardekooper & Papaloizou 2009a). 

It is still quite uncertain how the horseshoe drag de- 
saturation operates in radiative disks. The vortensity- 
related horseshoe drag presumably requires the same con- 
straint on viscosity as in isother mal disk models. Simi- 
larly, Baruteau & M assed (l2008ah argued that the entropy- 
related horseshoe drag should remain fully unsaturated if 
the timescale for restoring the unperturbed entropy gradi- 
ent is also bound by the horseshoe U-turn and libration 
ti mescales. This result has been i ndepend ently verified 
by IPaardekooper & Pa paloizoul ( 120081) and by 'K iev & Cridal 
(120081) . who respectivelv considered thermal diffusion and ra- 



diative effects as the dominant physical process sustaining the 
entropy gradient. However, the desaturations of both parts 
of the horseshoe drag are likely to be coupled. In particular, 
some viscosity should always be required to unsaturate the 
entropy -related horseshoe drag, so as to allow the exchange of 
angular momentum between the horseshoe region and the rest 
of the disk. A systematic study involving different viscous 
and thermal diffusion coefficients will certainly give more in- 
sight into the desaturation process in radiative disks. 

The horseshoe drag also plays a key role in another mech- 
anism potentially halting migration. Masset et al. (200<)l) 
showed that type I migration can be stalled near positive sur- 
face density transitions, when planets migrate toward under- 
dense regions in the protoplanetary disk. These density tran- 
sitions can be caused by ionization transitions, which may oc- 
cur near the snow line ( Kretke & Lin 2007), or near the in- 
ner edge of the dead zone (the region near the disk mid-plane 
sandwiched together by partially ionized surface layers). As 
the planet drifts inwards, the slope of the local density profile 
becomes increasingly positive. The positive horseshoe drag 
strongly increases, and migration halts when it balances the 
negative differential Lindblad torque. The ability of these so- 
called planet traps to stall migration depends on the horseshoe 
drag amplitude, hence on the local dissipation processes reg- 
ulating its desaturation. 

1.3. Turbulence 

Whether one invokes the gas thermodynamics to slow 
down, halt or reverse migration, or the presence of planet 
traps to stall migration, dissipation processes, including vis- 
cosity, must be taken into account. In the framework of planet 
migration, protoplanetary disks are commonly assumed to be 
laminar, with a constant viscosity whose amplitude aims at 
modeling the turbulent transport of angular momentum. 

Turbulence may originate from hydrodynamic instabilities, 
such as the Rossby-wave instability (Lovelace et al. 19991 
the global baroclinic instability (Kl ahr & Bodenheimein 
120031) . or the Kelvin-Helmholtz instability, triggered by the 
vertical shear of the gas as dust settles in the mid-plane 
(Johansen et al. 2006). Another source of turbule nce is the so- 
called magnetorotational instability (MRI, B albus & Hawleyl 
119911) . This instability relies on the coupling of the ionized 
gas to the weak magnetic field of the disk. Ionization may 
occur in the vicinity of the central object due to the star irradi- 
ation, or further out in the disk layers, most probably through 
the UV background or cosmic rays. The ionization state of 
the disk mid-plane, where protoplanets build up and migrate, 
is rather unclear It is still a matter of debate whether the disk 
mid-plane should be ionized as well, or if it remains neutral, 
being shielded by the layers. In an y case, its ionization st ate 
may significantly evolve with time ftlgner & Nelsonll200 8). 

So far, only magnetohydrodynamic (MHD) turbulence has 
been extensively studied in numerical simulations of planet- 
disk in teract ions . Nelson & Papaloizou (200 3). Wint ers et al.l 
(120031) and iPapaloizou et al. (2004j) have investigated the 
properties of gap form ation b y gia nt planets in tur bulent disks. 
iNelson & P apaloizou (12004 and iNelsonI (120051) have revis- 
ited type I migration for one or several planets embedded in 
turbulent disks. They found that low-mass planets evolve on 
a random walk, in contrast to the systematic decay expected 
in laminar disks. It is unclear whether, with longer simula- 
tions, the time-averaged stochastic torque due to turbulence 
would become negligible compared to the total mean torque. 
Another open question is whether this total mean torque sig- 
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nificantly differs from that obtained in laminar disks. In this 
regard, a systematic study of the horseshoe drag behavior with 
turbulence has not been done yet. We eventually point out that 
all aforementioned studies have assumed a fully magnetized 
disk. Simulations of planets interacting with disks harboring 
a dead zone have not been performed yet. 

1.4. Paper outline 

In this paper, we investigate the impact of turbulence on the 
horseshoe drag. Can the horseshoe drag remain unsaturared 
with turbulence? Can migration be stalled at planet traps with 
turbulence? Three-dimensional MHD calculations would cer- 
tainly be an elegant way to answer these questions. Nonethe- 
less, the horseshoe drag saturation typically occurs on several 
hundreds of orbits for planets subject to type I migration. The 
computational cost of 3D MHD calculations properly resolv- 
ing the horseshoe region over several hundreds of orbits is 
prohibitive. We perform instead two- dimensional hydrody - 
namic simulations using the model of iLaughlin et al.l (l2004t) 
(hereafter LS A04) for mimicking the properties of MHD tur- 
bulence. Further simplification is obtained by assuming an 
isothermal disk. The inclusion of the disk thermodynamics 
will be investigated in a subsequent study. 

The plan of the paper is as follows. In §|2] we describe the 
numerical set-up of our simulations. The turbulence model 
used in our study is detailed in § |3] in particular we assess the 
minimum duration of the calculations to reach torque conver- 
gence. Our calculation results on the horseshoe drag desatu- 
ration with turbulence are presented in § |4] Discussion and 
conclusions follow in § |5] 

2. NUMERICAL SET-UP AND NOTATIONS 

To investigate the desaturation properties of the horseshoe 
drag with turbulence, we perform two-dimensional hydrody- 
namic simulations using the turbulence model of LSA04. Be- 
fore detailing the physical properties of this potential in § [3] 
we describe in this section the numerical set-up of our calcu- 
lations. 

We use the two-dimensional code FARGO^^ in its isother- 
mal version. It is a staggered mesh code that solves the 
Navier-Stokes and continuity equations on a polar grid. An 
upwind transport scheme is used along with a harmonic, 
second-order slope limiter (Ivan Leeiill977l) . Its specificity is 
to use a change of rotating frame on each ring of the grid, 
which increases the timestep significantly (Masset 2000a b). 

The code units are the following. The initial orbital radius 
Vp of the planet is the length unit, the mass of the central object 

is the mass unit, and (GM*/rp^)^'/^ is the time unit, G 
being the gravitational constant (G = 1 in our unit system). 
The orbital period at r = is denoted by Toib. 

We denote the polar coordinates by r and q). The disk is 
initially axisymmetric and rotates at the sub-Keplerian an- 
gular velocity i2(r). The disk is isothermal, its aspect ra- 
tio reads h{r) = hp x (rjrpf-l'^, with hp ranging from 3% 
to 7%. The initial surface density is Ep x (rjrp^^^, with 
Ep = 5 X lO^'* and a = 0.5 by default. The smallest value 
of the Toomre parameter at r = rp then amounts to Qp ~ 20. 
The disk self-gravity is thus neglected throughout this study. 
Should we have included the full self-gravity, the differen- 
tial Lin dblad torque would have been larger by less than 5% 
(iBarutea u & Masset 2008b, figure 8). No kinematic viscosity 

3 See: |http://fargo.in2p3.fr | 



is included when the turbulent potential is applied to the disk. 
For comparison, runs without turbulence, but with a constant 
kinematic viscosity v are also performed. They will be re- 
ferred to as laminar runs, whereas calculations featuring the 
turbulent potential will be mentioned as turbulent runs for fu- 
ture reference. 

The planet mass is denoted by Mp, and q is the planet to 
primary mass ratio. The planet's Bondi radius rs is defined 
as rg ~ GMp/c1, where is the sound speed at the planet 
location. The softening length of the planet's potential is 
e = 0.6hprp, so that rg/e w 1.7(^//ip). For the planet masses 
considered in §|4l r^/e ranges from 0.3 to 0.6. 

The grid used in our calculations has A^^ = 5 12 radial zones, 
and Ns — 1536 azimuthal zones. The disk extends from 
'"min — 0.4 Tp to Tmax = l-8rp along the radial direction, and 
wave-killing zones are used next to the boundaries to mi ni- 
mize unphysical wave reflexions (Ide Val-Borro et al]|2006h . 

3. TURBULENCE PROPERTIES 

The turbulence model used in our study is that of LSA04. 
It is based on applying an external turbulent potential to the 
disk. We show in this section that the properties of this po- 
tential (spatial and temporal fluctuations, amplitude) can be 
tuned so that the perturbations it induces much resemble those 
obtained with 3D MHD calculations. As in the original work 
of LS A04, the turbulent potential 't'turb applied to the disk cor- 
responds to the superposition of 50 simultaneous wave-like 
modes. It reads ( lOgihara et 10120071 hereafter OIM07): 

50 

^mrh{r,^,t)^yr^D}Y^Ak{r,^,t), (1) 

k=l 

where 7 is a dimensionless constant indicating the turbulence 
strength, and 

= ^i^e"*'""''*'"/'^* cos{m{k)(p -(p^- Q-ktu) sm{ml/Atk). 

(2) 

In Equation (|2]i, {rk:^k) denotes the initial location of the 
mode of wavenumber m{k). Both r^^ and <pi^ are randomly 
sorted with a uniform distribution. The modes radial extent 
is Ok — nrii/4m. Modes start at time to j^, and their lifetime is 
Atk — iTZr^/mCs, with Cs the local sound speed. We denote by 
Q.k the Keplerian angular frequency at r = r^, fi^ = t — Jq and 
is a dimensionless constant sorted randomly with a Gaus- 
sian distribution of unit width. Note that the expression of the 
turbulent potential is independent of the disk surface density. 

3.1. Power spectrum and modes truncation 

LSA04 randomly sorted m with a logarithmic distribution 
between m = 2 and m ~ A^sec/8. Following OIM07, we 
also include m = 1 modes, and we set A<. = if m{k) > 6. 
The latter assumption is motivated by the fact that high-m 
modes contribute less to the turbulent potential, since Aj^ oc 
exp(— m^), within smaller time intervals (Af^. 0= 1/m). They 
therefore contribute less to the turbulent torque exerted on the 
planet. As we shall see hereafter, excluding m > 6 modes 
helped reduce the convergence time of our calculations. This 
cut-off impacts the power spectrum of the density perturba- 
tions triggered by the turbulent potential. The amplitude c„, 
of the corresponding Fourier coefficients, which read 

_ \j:-Ci:ir,cpj)rdrdcpe-"-1>\ 

T'liiax T-^'^V/,. +\ 1 ^ 
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Fig. 1. — Time-averaged power spectrum of the density perturbations due 
to turbulence. Results are shown with the m > 6 modes truncation (stars), 
and without it (diamonds). Both calculations have 7 = 4 x 10^^. No planet 
is included. 



Fig. 2. — Autocorrelation function of the torque per unit mass exerted by 
the turbulent disk. We compare the autocorrelation timescales obtained with 
the modes standard lifetimes (LSA04, stars) and reduced lifetimes (standard 
lifetimes divided by ten, diamonds). 



were calculated for two simulations with 7 = 4 x 10^^, and 
with no planet. One includes the m > 6 modes truncation, 
the other does not. We display in Figure [T] the coefficients 
Cm time-averaged over the runs duration (about 700 Torb). The 
density spectrum without truncation, and that of the 3D MHD 
simulation of LSA04 (see their figure 2), have very similar 
slopes. The upper dashed line overplotted in Figure [T] indi- 
cates that this spectrum decreases as m^^/^. We comment 
that in absence of modes cut-off in the turbulent potential, 
the power spectrum features both the modes that are directly 
forced by the turbulent potential, and the modes induced by 
energy cascade from larger scales. A detailed comparison 
with Kolmogorov theory of weakly compressible fluids tur- 
bulence would not be strictly relevant here. Spectra with and 
without the m> 6 modes cut-off much resemble up to m = 6, 
below which the spectrum with truncation drops off, and 
starts decreasing as m^^. This is the slop e expected in two - 
dimensional forward decaying turbulence dKraichnanll 1 9671) . 
The turning point from m w 30 is not expected, however. With 
increasing the turbulence amplitude to 7= 10^^, we find that 
the turning point moves to a higher value (m w 60), while the 
spectrum without truncation is not significantly altered. Ad- 
ditionally, we checked that the spectrum of velocity perturba- 
tions is very similar to the density spectrum, as expected with 
an disk. Unless otherwise stated, all calculation results shown 
below are obtained with the m> 6 modes cut-off. 

3.2. Autocorrelation timescale 

The autocorrelation timescale Tc measures the typical 
timescale over which perturbat ions due to tu rbulence become 
uncorrected. Following Ois hi et al.l (l2007h . we evaluate the 
autocorrelation function of the torque per unit mass, F, ex- 
erted by the turbulent disk on a massless planet. Its normal- 
ized value is given by 



ACF(t) 



S'r^T{t)T{t-x)dt 

j^"'--"'r^{t)dt ' 



(4) 



where fmax is the runs duration. The quantity ACF is displayed 
in Figure |2]for two calculations, one with the modes lifetimes 
taken by LSA04 (Afj. ~ Inrj^/mcs, they are referred to as stan- 
dard lifetimes), and one for which Af^ is reduced by a factor of 



ten (reduces lifetimes). For both runs, 7 = 4 x 10^^, and the 
disk aspect ratio at r = rp (planet's location if any) is hp = 3%. 
As in all other calculations, the sampling time of the torque is 

rorb/20. 

There are several methods for evaluating T^. One 
is based on calculating Tc as tc — /q™" ACF(t)c/t (e.g. 
iRein & Papaloizoul 12009*). With the standard lifetimes (de- 
noted by stars in Figure |2|, we find « 2roib, whereas 
Tc ~ 0.3 Torb with the reduced lifetimes (diamonds in the same 
figure). Alternatively, t^- can be estimated as the smallest lag 
value for which the autocorrelation function cross es zero with 
a positive slope (second zero-crossing, see e.g. lOishi et aP 
(2007)). From Figure |2] this method yields Tc « 6 Torb 
with the standard lifetimes, and w 0.8 Torb with the re- 
duced lifetimes. The autocorrelation timescale obtained with 
the reduced lifetimes is in good agreement with the typi- 
cal values obtained with 3D MHD calculations, be they lo- 
cal (Fromang & Papaloizou 2006; Oishi et al. 2007) or global 
jNelson & Pat)aloizoin i2004t jFromang & Nelson 2009), and 
be they obtained with disks fully i nvade d by the MRI, or in- 
cluding a dead zone (lOishi et al.l l2007h . Unless otherwise 
stated, the runs presented hereafter use the modes reduced 
lifetimes: Af^. = Q.lnrj^/mcs. 

3.3. Reynolds stress parameter 

Another relevant property of the turbulent potential is the 
transport of angular momentum it generates. The turbulence 
amplitude 7 can thus be related to the magnitude of a in 
the widely a dopted alpha prescription for turbulent viscosity 
dShakura & Sunvaev 1973). A usual procedure is to calculate 
the Reynolds stress parameter aR{r,t), which we evaluate as 



aR{r,t) 



(5) 



where £ denotes the axisymmetric surface density, P = Ec^ 
is the gas pressure, 5vv = v,- — v7 and 5v(p = — V^, with 
Vr (V(j)) the radial (azimuthal) component of the gas velocity. 
Here the overbar symbols denote surface density-weighted az- 
imuthal averages. The numerator in Equation (|5]) features the 
Reynolds stress tensor We furthermore average the profile 
aR{r,t) (i) with radius over the whole disk, then (ii) with time 
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Fig. 3. — Average Reynolds stress parameter («/;) calculated for several 
values of the turbulent coefficient y. For all runs, the disk aspect ratio at 
r = is hp = 3%. The solid curve depicts a quadratic function of / passing 
through the point (7=4x 10""*,(aR) = 6.2 x 10"^). 

over the calculations duration. The resulting number is de- 
noted by {(Xr). We performed a series of runs with varying 7, 
following the set-up described in §|2](/ip = 3%). We display 
(a«) as a function of 7 in Figure[3] The solid curve highlights 
that (ur) scales as 7^. This scaling is expected as both the 
perturbed velocities 5v,- and 5v(p scale with 7. 

We investigated the dependence of (ur) with various pa- 
rameters of the turbulent potential. Firstly, we found that, for 
a fixed value of 7, (ur) is typically increased by a factor ^ 1 .5 
without the m > 6 modes truncation. This confirms that most 
of the angular momentum exchange driven by turbulence oc- 
curs through (very) low wavenumbers. Secondly, we varied 
Tc by setting the modes lifetimes to Ati^ = nx Inri^/mcs, with 
n varying from 0.02 to 2. By taking the second zero-crossing 
of the torque autocorrelation function, we measured in the 
range [0.6 — 8] Tdb. Note that does not scale with n. We 
found that (a^) scales with for Xc < TLrb. Beyond, {ur) 
slightly decreases with Tc (by ~ 20% when increasing from 
w 2rorb to w STorb). Lastly, another series of runs with sev- 
eral values of 7 was performed for hp — 7%. Altering hp (or, 
equivalently, the sound speed) modifies the pressure (P /i^) 
and the modes lifetimes (Atk °= hj,^). From the comparison 
with the previous series for hp = 3%, we found that (afl) 
scales with /Zp^, suggesting the decrease of the modes life- 
times was here of negligible importance. We checked indeed 
that reducing hp from 7% to 3% led to decrease by only a 
~ 20% factor. The results obtained with both series of runs 
therefore lead to the following relationship: 

(afi)«35(^^^ , (6) 

which inverts as 7 w 1 .7 x 10"^ (ocr). 

3.4. Diffusion coefficient 

The saturation level of the horseshoe drag depends on the 
advection-diffusion of vortensity inside of the planet's horse- 
shoe region. In laminar disks, vortensit y dif fusion is con- 
trolled by the kinematic viscosity (see § ll.2l i. while in tur- 
bulent disks it should depend on the turbulence strength. In 
§ 13.31 we have related the turbulence amplitude 7 to an equiv- 



FlG. 4. — Vortensity profiles obtained for laminar and turbulent runs with 
the same equivalent Reynolds alpha parameter («/;). Both profiles are sub- 
tracted from the background profile. A slight negative perturbation (10% 
relative difference) is initially imposed to the unperturbed vortensity pro- 
file at r = r/j, the location of which is depicted by the vertical dotted line. 
Stars show the perturbed vortensity profile at t = 4.5 To^b for a laminar run 
with V = 3.4 X 10"'. The dash-dotted curve shows the profile expected from 
Equation for D = v. Diamonds display the profile averaged over four tur- 
bulent calculations with 7= 10"*, and time-averaged between 4 and 5 orbits. 
The solid curve shows the result of Equation 0, also time-averaged between 
4 and 5 orbits, with D = 4v. 

alent Reynolds alpha parameter (Ur), which quantifies the 
turbulent transport of angular momentum. There is no rea- 
son to expect that the turbulent viscosity associated with (Ur), 
namely the quantity Vr = {aR)csH, should correspond to the 
vortensity's diffusion coefficient, which we denote by D. To 
evaluate D, we altered the unperturbed surface density profile 
by imposing a slight overdensity (10%) on one ring, located 
at r = Til. A perturbation of same amplitude, but of opposite 
sign, is thus applied to the unperturbed vortensity profile. In 
absence of viscosity and turbulence, the density profile flat- 
tens out while the vortensity profile is maintained at its initial 
value, to within the effects of numerical viscosity. This actu- 
ally allowed us to check that the numerical viscosity is much 
smaller than the minimum values of laminar and turbulent vis- 
cosities used in our study. For laminar and turbulent runs, 
we checked that the time evolution of the vortensity profile V 
(subtracted from the background profile) takes the form 

where Vd is the initial vortensity perturbation at r = r^, and 
5r is the mesh radial resolution. This is illustrated in Fig- 
ure |4] for a laminar run with v = 3.4 x 10"^, and a turbulent 
run with 7= lO"'*. For both calculations, hp = 3%, ~ 1.1, 
and IVd'Lp/ Q,p = 0.1. The values of v and 7 were taken such 
that (aff), given by Equation is equal to v/csH. Both 
calculations thus have the same equivalent Reynolds alpha 
viscosity: (ur) ~ 3.8 x 10"^. The vortensity profile (sub- 
tracted from the background profile) is shown at t = 4.5 Torb 
for the laminar run (stars). The dash-dotted curve is the per- 
turbed profile expected from Equation ^ with D = v, at the 
same time. It shows that, as expected, the vortensity's dif- 
fusion coefficient i n laminar disks correspo nds to the kine- 
matic viscosity (e.g. lOgilvie & Lubowll2006l) . Also displayed 
is the perturbed vortensity profile time-averaged between 4 
and 5 orbits, and further averaged over four different turbu- 
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lent calculations (diamonds, the turbulent runs differ by the 
random numbers sorted in the expression for the turbulent po- 
tential). The solid curve shows the perturbed profile given by 
Equation (|7]i, time-averaged over the above time range, with 
D = 4{ai{)csH = 4v. The good agreement with the numerical 
profile underscores that, in our turbulent tuns, the vortensity's 
diffusion coefficient differs from the turbulent viscosity Vr as- 
sociated with the turbulent transport of angular momentum. 
By considering several other values of 7, we have estimated 
the ensemble average of the vortensity's diffusion coefficient 
as D w 4-Vr, with a relative uncertainty of w 30%. Differently 
stated, vortensity diffusion is four times more efficient than 
angular momentum transport in our turbulent calculations. 

The ratio Vr/Dis usually referred to as the Schmidt number 
Sc- For the vortensity radial diffusion, we thus find a (radial) 
Schmidt number of Sc ~ 0.25. Similar values, albeit slightly 
smaller, were obtained with the diffusion of passive scalars. 
For comparison, studies with 3D MHD calculations report ra- 
dial and v ertical Schmidt numbers f or dust diffusion of order 
unity (seel Fromang & Nelsonll2009 , and references therein). 
Note that in such simulations, Vr is substituted by Vr + Vm, 
where Vm is the viscosity associated with the Maxwell stress 
tensor. Schmidt numbers are however sensitive to the parti- 
cles sizes, and they can become much smaller than unity for 
small p articles well-coupled to th e gas, thus acting as passive 
scalars dPromang & Nelsonll2009l) . 

In laminar disks, the kinematic viscosity v is commonly 
modeled by a dimensionless alpha viscosity a = v/csH. Sim- 
ilarly, we define for turbulent runs an equivalent alpha vis- 
cosity (ao) as (od) = D/csH. Note that {ao) = S^^ccr) w 
4{aR). Equation (|6| can thus be recast as 



Y ~ 8.5x10-^ hp ^(od) 



(8) 



which inverts as (ud) ~ 1.4 x lQ^{Y/hp)^. Interestingly, in 
MHD calculations, the Maxwell stress parameter is typi- 
cally threetimesj^;gerthan the Reynolds stress parameter Ur 
(e.g. lFromang & Nelsonl [2009). The total alpha parameter a 
thus verifies a ^ 4-aR. Since the Schmidt number is also com- 
monly of order unity in MHD calculations (but see above), 
we suggest that our value of (ao), given by Equation (|8]l, can 
model the total alpha parameter of MHD simulations. 

Typical alpha parameters in active zones are in between 
10^^ and 10^\ depending on the geometry of the magnetic 
field (see e.g. Fromang & Papaloizou 2007). Even if the mag- 
netorotational instability only develops in the disk laye rs, tur- 
bulen ce may be induced in the dead zone (Fleming & Stone| 
l2003h . These authors estimated that the alpha parameter in 
dead zones might be about two orders of magnitude smaller 
than in active layers. Values of a ranging from 1 0~^ to 10^** 
might thus be relevant for dead zones (see also lOishi et aTl 
|2007). From Equation ([8]), and assuming hp — 5%, we find 
that values of 7 in the range [10^''^ — 10^^] should be relevant 
for active layers, while 7 < 5 x 10"^ should be more appro- 
priate for dead zones. 

3.5. Running-time average of the turbulent torque 

In addition to potentially desaturating the horseshoe drag, 
the turbulence also exerts a stochastic torque on the planet 
through the density perturbations it generates. This stochas- 
tic torque is referred to as the turbulent torque and is de- 
noted by Ftuit- We are primarily interested in the cumula- 
tive effects of the turbulent torque, which we can evaluate 
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through its running- time average Fturb, defined as FturblO = 
f^' X Jq rturb(M)i^M. To investigate the properties of Fturb, we 
performed a series of calculations with a massless planet, and 
with varying the turbulence coefficient 7. The specific turbu- 
lent torque (turbulent torque by unit planet mass) is directly 
given by the instantaneous specific torque exerted by the disk 
on the planet. We display in Figure |5] the distribution of the 
specific turbulent torque for 7 ~ 10^^, after 1000 orbital pe- 
riods. A Gaussian best fit is overplotted and shows that the 
turbulent torque frequency follows a Gaussian distribution. 
Thu s, the amplitude of Fturb can be written as (e.g. .Nelsoiil 
5005h 



I Fturb I ~ O'turb X 



-1/2 



for f > Tt, 



(9) 



with (Tturb the mean deviation of the turbulent torque distri- 
bution. It takes the form CTtmb = CUpgyrpSlj,, where C is 
a dimensionless constant, and where the quantities with a p 
subscript are to be evaluated at the planet location (LSA04, 
OIM07). It can be recast as CTturb = ClnGLpMprp, which 
corresponds to a fraction C of the natur al scale for the torque 
exerted by the disk on the planet (e.g. iJohnson et alJl2006) . 
The functional dependence of (Tturb is expected from the ex- 
pression of the turbulent potential, given by Equation (|T|). In 
particular, aiurb does not d epend on the disk scale height, as 
in 3D MHD calculations (Nelson & Papaloizou 2004j). 

We checked with our series of calculations that aturb scales 
with 7. The value of the constant C was inferred from the 
turbulent torque distributions. We find that: 



CJturb«2.4x 10"Ep<77r;n; 



(10) 



when the m > 6 modes truncation is taken into account. The 
above expression should otherwise be multiplied by a factor 
of w 1.25. This is about 20 times larger than the expression 
given by LS A04. We show in § 14.41 that such large differ- 
ence cannot be explained by the absence of m = 1 modes in 
the work of LSA04. It is also interesting to compare the ex- 
pression in Equation ( fTOl i with the results of 3D MHD cal- 
culations. Using Equation ([8]) together with typical disk pa- 
rameters in the 3D MHD simulations of .Nelson (.2005.) - 
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a w 8 X 10-^ at r,, « 2.5, hp = 7%, and Upr,, w 5.2 x IQ-^^ 
jJohnson et al.ll2006l) - we find Oturb/q ~ 7 x 10^^, which is 
not far from the value found by "Nelson (2005) 1 .5 x 10""* 
in the same code units, see his figure 14). This general agree- 
ment makes us confident that the results obtained with our 
two-dimensional simulations reproduce the main turbulence 
properties of realistic 3D MHD calculations. We finally com- 
ment that Equations (|9]l and ( fTOl i have been worked out for a 
massless planet. The turbulent torque actually felt by a low- 
mass planet may be different, as the disk responses to the 
planet potential and to the turbulent potential could be cou- 
pled. 

3.6. Convergence time 

We assume that the torque on the planet is the sum of the 
differential Lindblad torque AFlr, the horseshoe drag Fhs, 
and the turbulent torque Fturb- The results shown in §|4]are 
notably aimed at testing this assumption. To properly assess 
the desaturation level of the horseshoe drag with turbulence, 
we must ensure that | Fturb | is a small fraction of | AFlr + Fhs |- 
We denote by / this fraction. This constraint provides the 
convergence timescale fconv of the calculations. From Equa- 
tion (|9]l, we have 

W^T,/-^xf (11) 
VAiLR + iHS/ 

The differential Lindblad torq ue is ev aluated by a recent es- 
timate derived by Pa ardekooper et al.l (12009). This estimate, 
which updates the standard formula of lTanaka et al.l (120021) by 
including the softening length of the planet potential, reads: 

AFLR = -CLR?2Epr4n2/Zp^ (12) 

with 

Clr«(2.5 + 1.7t-0.1(7)x f— -^j , (13) 

where T = — t/logro/c/logr is the power-law index of the 
unperturbed temperature profile Tq, and Hp = h„rp. Equa- 
tion (fT2] | differs from the 3D estimate of lTanaka et al.l (120021) 
by the constant Clr, which reads Clr « 2.340 - 0.099a + 
0.418t in the latter work'*. For isothermal disks (t — 0), both 
formula give very similar results, as far as conservative values 
of the softening length are concerned. For £ /Hp = 0.6, they 
differby-25%. 

For the sake of simplicity, we assume that the horseshoe 
drag takes its fully unsaturated value. I n isothermal disks, the 
fully unsaturated horseshoe drag reads (IWardlll99l1) : 

riis = l(^l-(yyt^pal, (14) 

where is the half-width of the planet's horseshoe region, 
and where 3/2 — (7 is the opposite of the vortensity gradient. 
By equating the expression in Equation ( fT4] i with the two- 
dimensional linear corotation torque formula of lTanaka et alj 

We comment that lTanaka et all 420021) provided estimates of the differen- 
tial Lindblad torque, and of the linear corotation torque for locally isothermal 
disks (except for the linear corotation torque in two dimensions). In three 
dimensions, the total linear torque reads Fsd = —C3u<]^^p''pClj,hp-^, with 
C3D = 1.364 + 0.541C7 + 0.433T. 



( 120021) . iMasset et all ( l2006ah estimated the half-width of the 
horseshoe region as 

x,^l.lrp[^y'\ (15) 

Usin g a simple tw o-dimensional model for the horseshoe re- 
gion, |P^d^rooper & Papaloizou ( 2009b) have recently de- 
rived an analytic expression for featuring the softening 
length of the planet's potential (see their equation 39). For 
small softening (e ^ O.lHp), their expression agrees with 
the values of Xs measured in numeric al simulation s , and i t 
differs from the above expression of IMasset et all (l2006ah 
by a factor of ~ 1.5, quite independently of e (see their 
figure 10). For e > 0.1 //n, the analytic expression of 
iPaardekoop er & Pap aloizoul (^009b) is reduced by a fitting 
factor of ^ 30%, meant to reproduce results of simulations. 
The resulting estim ate is in good agreement with that of 
IMasset et al.r(l2006al) for e/Hp ranging from 0.2 to 0.6. For 
the run parame ters of jj|4l(<7 //i^ ^ 0.2, e — 0.6Hp), they agree 
to within 10% ( Paardekooper & Papaloizoul20 09b, figures 10 
and 11). They are also in good agreement with our results 
of simulations, which we checked by a streamline analysis. 
Using Equation ( fTSl l, Equation ( fT4] i can then be recast as 

Tiis = -CHsq^^p4^lhf, (16) 

with 

Chs«-1.1x (17) 

Combining Equations (O, ( fTOj ) to ( fT3] l, ( fT6] l and (fTTT l, we 
are left with: 

where we recall that the horseshoe drag Fhs was assumed to 
take its fully unsaturated value. Equation ( fTSI l may be refined 
by including the analytic dependence of Fhs with viscosity 
expected in laminar disks (Masset 2002). This would come 
to substituting Chs in Equation (VH by Chs -^{zs), where the 
quantity J^(zj) is given by Equation ( |22] |. When Chs < 0, 
namely for a < 3/2, the expression in Equation (fTSl l corre- 
sponds to a maximum convergence time. 

Further insight can be obtained by estimating the range 
values of (ao) such that the horseshoe drag is unsaturated. 
As recalled in § 11.2! this condition is roughly fulfilled in 
laminar disks when the diffusion timescale across the horse- 
shoe region, Tyisc ~ x^/v, is smaller than the libration period 
Tiib = ^Kr„/ (3Q.oXs), but larger than the h orseshoe U-turn 
time Tu_tum 

« /jpTiib (Baruteau & Masset 2008a). Anticipat- 
ing that this condition essentially holds with turbulence (see 
I4.3l l. and assuming that Xs ~ l.lrp{q/hpYl^, we find that 
the value of {ao) air ~ rp should verify 

Q.l6q^l^hp'^'^ < (ao) < Q.ieq^l^hp'^'^. (19) 

We therefore assume that the horseshoe drag is maintained 
at its fully unsaturated value for {ud) ~ Q.\6q^l^ h/j'^ . This 
estimate is checked against numerical simulations in § 14.21 
Equation ( flSl ) can finally be recast as 

fconv ~ 65 T, — ^ ^ . (20) 

/2 Clr + ChS ' 
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We emphasize again that the expression given in Equa- 
tion ( l20b assumes that the horseshoe drag is fully unsaturated, 
and that Chs is independent of e. For our purposes, minimiz- 
ing the convergence time implies a compromise between a 
small aspect ratio, and a planet mass that is not too large to be 
relevant for type 1 migration. The disk andplanet parameters 
used in our calculations are described in § |4] 

4. RESULTS OF NUMERICAL SIMULATIONS 

We present in this section the results of simulations with 
planet, and with the turbulence model detailed in § [3] These 
simulations follow the numerical set-up described in § |2] Re- 
call in particular that the power-law index of the initial sur- 
face density profile is a = 0.5. The disk and planet param- 
eters were chosen to minimize the maximum convergence 
time Tconv of the turbulent runs, given by Equation ( |20] i. We 
took hp — 3% and ^ = 5 x 10^^, which corresponds to a 
Mp « I.6M0 planet mass if the central object has a solar 
mass. Using Tc ~ 0.5 orbits (see ^ \3.2\i . and assuming / = 0.1, 
we have Tconv ^ 2000 orbits. Our simulations were run for 
about 4000 orbits, so we expect that iFturbl should not exceed 
~ 0.1|AFlr + Fhs| at the end of the simulations. For com- 
parison, laminar runs without turbulence, but with a constant 
kinematic viscosity V are performed. Runs with turbulence 
required about twice as much computing time as runs with- 
out turbulence. We also point out that the half-width of the 
planet's horseshoe region is ~ 0.015rp. It is resolved by 
about 6 cells along the radial direction. This resolution is sim- 
ilar to that of other recent numerical studies of the horseshoe 
drag. The i nflue nce of the grid resolution on our results is 
assessed in § 14.41 

4. 1 . Torque vs. turbulence amplitude 

We performed a series of 9 runs with varying the turbu- 
lence amplitude from 7= 10^^ to 7= 1.5 x 10^^. Using 
Equation (O, the equivalent alpha viscosity (ao) related to 
the vortensity's turbulent diffusion ranges from 1.5 x 10^^ to 
3.5 X 10^^. A laminar inviscid run (7 = 0, V = 0) was per- 
formed for comparison. The top panel of Figure|6]displays the 
relative perturbation of the disk's surface density (with respect 
to the initial profile) obtained at 500 orbits for the laminar in- 
viscid run. A shallow gap is progressively opened up, which 
is consistent with the non-linear wave dissipation model of 
iRafikov (2002) in inviscid disks. This model predicts that the 
minimum planet to primary mass ratio ^niin for gap opening is 

?min = ^ X Min{5.2ep^/\3.8(e„Ap)"''"'} • (21) 

The Toomre parameter at the planet location is Qp ^ 20, and 
Equation (l2Tl l gives w 5.7 x 10^^. This value is in agree- 
ment with our planet to primary mass ratio. Note however that 
this model does not include the planet's softening length, and 
the possible interaction betwe en the planet w ake and the flow 
inside of the horseshoe region ( iMasset et al.l2006 a). We men- 
tion that a shallow gap also formed for the lowest values of the 
turbulence amplitude. The gap impact on the torque variation 
remains small, as will be shown below. The middle panel of 
Figure |6] displays the perturbed density for the turbulent run 
with Y — lO^'*, which corresponds to {ao) ~ 1.5 x 10^^. For 
comparison, the perturbed density obtained with 7 = 10^^, 
but without the m> 6 modes truncation, is shown in the bot- 
tom panel of Figure |6] In both cases, the planet wake and 
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Fig. 7. — Time evolution of the running-time averaged specific torque exerted on a Mp « \.6Mq planet mass. Results ai'e shown for several values of y, for 
a = 0.5 (left panel) and fj = 1.5 (right panel). The dashed and dash-dotted lines in the left panel depict the values of the differential Lindblad torque and of the 
fully unsaturated torque, respectively, both measured in a run without turbulence. 



the density perturbations due to turbulence have similar am- 
plitudes, and no gap is visible. For these runs, the averaged 
Mach number associated with the velocity's turbulent pertur- 
bations is w 0.3. We comment that a similar wave pattern 
is obtained in absence of any embedded planet. Note that 
the spiral waves generated by our turbulence model are more 
tightly wound close to the disk's inner edge, since the modes 
radial extent scales with radius. 

The time history of the running-time averaged torques 
(hereafter, r.t.a. torques) obtained with our series of runs are 
depicted in the left panel of Figure |7] The torques seem to be 
all converged with time by the end of the simulations. The 
dashed and dash-dotted lines show respectively the values of 
the fully saturated and unsaturated torques, both measured in 
the run without turbulence. We comment that the fully satu- 
rated torque experiences a slow, stationary increase due to the 
continuous building up of a gap around the planet. For our set 
of planet and disk parameters, the fully saturated torque only 
increases by ^ 4% between 500 and 4000 orbits. It can thus 
be confounded with the differential Lindblad torque. As / in- 
creases, the stationary value of the r.t.a. torque increases from 
the differential Lindblad torque up to the fully unsaturated 
torque (here for 7 w 8 x 10^^), before it decreases with 7. 
This behavior reminds us of the torque dependence with vis- 
cosity in laminar disks (lMassetll2002h . This result is a strong 
indication that turbulence can unsaturate the horseshoe drag, 
depending on the turbulence strength. 

To further investigate this result, we performed another se- 
ries of calculations with an initial surface density profile de- 
creasing as r^^l^. In laminar isothermal disks, the horseshoe 
drag cancels out at all time, and the total torque is equal to 
the differential Lindblad torque. The results of these simula- 
tions are displayed in the right panel of Figure |7] for 7=0, 
7 = 2 X 10^^ and y — 10^^. Running-time averaged torques 
are depicted, even for the run without turbulence. For this 
run, recall that the slow increase of the torque is due to 
the gap clearance. The convergence time is clearly much 
smaller than in the previous series of runs with a = 0.5. Us- 
ing Equations (O and ( fTSl ). and assuming / = 0.1, we have 
fconv ~ SOOTorb (7/10^"^)^. For the values of 7 considered 
here, iFtuitl should be already smaller than ~ 0.1 IAFlrI from 
^ 500 orbits, which is clearly seen in the right panel of Fig- 



ure I2] The torques with turbulence rapidly converge toward 
the torque without turbulence. This result could be somewhat 
surprising since for the largest values of 7 of our study, the 
density perturbations due to turbulence are typically as strong 
as those triggered by theplanet wake (see the middle and 
bottom panels of Figure |6]l. It indicates that, for our range 
values of 7, and for the duration of our runs, the differential 
Lindblad torque is not significantly altered by turbulence, and 
that the torque on the planet can thus be decomposed into a 
stationary component (the differential Lindblad torque) and a 
fast-varying component (the turbulent stochastic torque), with 
no significant coupling between both. In addition, it justi- 
fies that the torque variation with 7 obtained with a = 0.5 
does arise from the horseshoe drag desaturation through tur- 
bulence. More insight into the l ong- term evolution of runs at 
high-turbulence is provided in ij 14.21 

4.2. Comparison to laminar runs 

We have shown in 14. II that turbulence can unsaturate the 
horseshoe drag. We now study how the torque dependence 
with the turbulence amplitude compares to that of the torque 
with viscosity in laminar disks. For this purpose, we per- 
formed an additional series of laminar runs with a = 0.5 
and a uniform kinematic viscosity v, which can be related 
to a dimensionless alpha viscosity a = v/csH. Recall that 
the vortensity's diffusion coefficient D in turbulent runs is 
similarly modeled by an equivalent alpha viscosity («/)) = 
D/csH. We shall assume hereafter that a and (ao) take their 
value at the planet location. For laminar runs, a varies from 
LI X 10^^ to 8.9 X 10^^. A steady state was attained between 
300 to 500 orbits. Stationary torque values obtained with the 
laminar runs are plotted against a in the left panel of Figure[8] 
Stationary running-time averaged torque values obtained with 
the turbulent runs are also depicted against 7 (top axis) and 
(«£)) (bottom axis, same scale as for a). 

The results of turbulent and laminar runs are globally in 
good agreement, especially for intermediate viscosity values. 
Slight differences at low viscosities are presumably due to the 
fact that torques were evaluated at different times (from 300 to 
500 orbits for laminar runs, between 3500 and 4000 orbits for 
turbulent runs). For the smallest viscosities, the r.t.a. torques 
are thus biased toward more positive values, triggered by the 
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Fig. 8. — Comparison of the turbulent and laminar runs for two planet to primary mass ratios: = 5 x 10 (left panel, turbulent torques are averaged over 
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depict the values of the differential Lindblad torque and of the fully u nsatu rated torque for an inviscid laminar ran. The sohd curve displays the dependence of 
the steady torque with a expected in laminar disks, given by Equation (22). 



progressive clearance of a gap around the planet. This bias de- 
creases with increasing viscosity. Note that the determination 
uncertainty in the relationship between 7 and («£>), given by 
Equation ([8]l, is another possible source of discrepancies. The 
overall good agreement between laminar and turbulent cal- 
culations at low-turbulence level implies that the total torque 
felt by a low-mass planet can be decomposed into a laminar 
torque and a stochastic torque due to turbulence. 

As in Figure I2I the dashed and dash-dotted lines show the 
values of the differential Lindblad torque AFlr, and of the 
fully unsaturated torque Fpu, both evaluated for 7=0 and 
V = 0. The solid curve depicts the analytic dependenc e of th e 
steady laminar torque with viscosity, given by Masset (i2001h . 
More precisely, we display the function T{zs) defined by 

Ffe) = AFlr + (Lfu - AFlr).^(z.v) , (22) 

where 

V g fc) / 

with Zs — Xs{2nah^i,r^p)^^ 1^ . In Equation ( l23T l. the function g 

is defined as g{z) — Bi(z) — V3Ai(z), where Ai and Bi de- 
note the Airy functions. The half-width of the horseshoe re- 
gion was determined numerically with the laminar run with 
a = 1.1 X 10^"* (dichotomic search of the separatrices, 
calculated as the geom etric average of t he ho rseshoe half- 
widths at ±1 rad, as in lCasoh & Massetl (l2009h '). We found 
Xs sa 1.429 X 10^^, which agrees to less than 1% with the 
estimate x,. w \.\rp{q/hpYl^ of iMasset elan (|2006a). The 
analytic expression of Equation (|22] | correctly reproduces the 
results of both the laminar and turbulent runs, for alpha vis- 
cosities smaller than ^ 0!,nax, where stationary torques take 
their maximum value. 

The differences between laminar and turbulent runs are 
more evident at high viscosities. Interestingly, laminar runs 
with a in the range [5 x 10^"* — 5 x 10^^] have steady torques 
that exceed the fully unsaturated torque obtained without vis- 
cosity. For any viscosity in this range, the torque excess ap- 
pears after a time comparable to the horseshoe U-turn time 



(which indicates that the torque excess is likely an excess of 
horseshoe drag), and it reaches a maximum value in about a 
libration time. The maximum laminar torque is obtained for 
Oimax ~ 2 X 10^^, which is in very good agreement with the es- 
timate {ao) ~ 0.16^^/^/1^7^ used to evaluate the convergence 
time at Equation ( |20l ). The relative difference between the 
maximum laminar torque, and the fully unsaturated inviscid 
torque is « 22%. Additional series of calculations^ revealed 
that this relative difference decreases with increasing aspect 
ratio: it amounts approximately to 10% for h = 5%, and to 
7% for h = 7%. A look at the horseshoe drag formulation of 
ICasoli & Massetl (120091) (their equation 22) suggests that the 
torque excess can arise from a modification of the horseshoe 
streamlines (e.g. a slight shift of the separatrices), or of the 
profiles of vortensity inverse T./(o at the upsti'eam parts of the 
horseshoe region (see e.g. Masset & Casoli 2009. their figure 
2). We performed a thorough streamline analysis that showed 
no significant modification of the horseshoe streamlines with 
varying viscosity. The profiles of E/o at the upstream parts 
of the horseshoe region are modified, however, in a way that 
we find to be consistent with the presence of a torque excess. 
The modification of these profiles with viscosity, and there- 
fore the mechanism that triggers the torque excess, deserves a 
detailed study, which goes beyond the scope of this paper. 

Turbulent runs do not feature such a torque excess. Torques 
with turbulence start decreasing at slightly smaller viscosities 
compared to laminar runs, and they seem to decrease faster 
To provide more insight into these differences, we show in 
the left panel of Figure |9] the surface density profiles, time- 
averaged over 4000 orbits, for some of the turbulent runs 
(solid curves). Density profiles obtained at 300 orbits for lam- 
inar runs with similar values of the vortensity's diffusion co- 
efficient are overplotted as dashed curves. At low-viscosity, 
laminar profiles are not strictly stationary, since a shallow 
gap ultimately forms around the planet, as can be seen in 
the turbulent run with 7 — 4 x 10^^. As already shown, the 

^ For these additional rans, we increased the resolution to keep the same 
ratio .«•,/ Sr, where Sr is the mesh size along the radial direction. 
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Fig. 9. — Surface density profiles (left panel) and vortensity profiles (right panel) for some of the turbulent and laminar runs displayed in the left panel of 
Figure[8] In both panels, profiles with turbulence are time-averaged over 4000 orbits (solid curves). Dashed curves show the profiles obtained at 300 orbits with 
laminar runs. For a fixed value of 7, laminar and turbulent runs have the same vortensity's diffusion coefficient. In the right panel, all vortensity profiles have 
been slightly offset to facilitate the comparison of their slopes. The vertical dash-dotted lines show the location of the separafiices of the horseshoe region. 



gap's impact on our results is weak, and we neglect it in the 
following. As 7 increases, profiles with turbulence tend to 
steepen around the planet, whereas laminar profiles almost 
coincide (and will stay so on the long term for the largest 
viscosities). Up to 7 = 10^^, laminar and turbulent profiles 
take similar values around the planet, which indicates that 
the time-averaged Lindblad torque with turbulence, AF^r, 
should approximately equal that of the laminar run, Ar^f,. 
For 7 = 1.5 X 10^^, the averaged density profile is reduced 
by w 15% with turbulence, which should decrease the Lind- 
blad torque accordingly, as the slope of the density profile has 
little impact on the Lindblad torque. To now estimate the im- 
pact of the density change with turbulence on the horseshoe 
drag, we show in the the right panel of Figure |9] the vorten- 
sity profiles for previous runs (time-averaged profiles for tur- 
bulent runs, and stationary profiles for laminar runs). They 
are depicted in a narrow range around the horseshoe region, 
the location of its separatrices are indicated as dash-dotted 
lines. Also, all profiles have been slightly offset to facilitate 
the comparison of their slopes. For 7 = 4 x 10^^, turbulent 
and laminar vortensity profiles hardly differ, which justifies 
that the total torques are in very good agreement. Nonethe- 
less, for higher values of 7, time-averaged vortensity gradi- 
ents are smaller with turbulence, by about 15%, 25%, and 
45% for 7 = 6 X 10"^ 7= lO""*, and 7= 1.5 x lO"'*, respec- 
tively. For 7 = 6 X 10^^, because the averaged density at the 
planet location is similar in turbulent and laminar runs, the av- 
eraged horseshoe drag with turbulence, F^g, should be sa 15% 
smaller than the stationary horseshoe drag of the laminar run, 
Fyg . Using the left panel of Figure [8] it turns out that, if 
the torque difference between laminar and turbulent runs was 
only due to the density change with turbulence, the averaged 
total torque with turbulence would be « AF^j^ + 0.85Fjjg = 
—3.06 X 10^^, which is about 7% larger than its actual value. 
Similarly, for 7 = lO^'*, the torque with turbulence would be 
w AF^R + 0.75F[js -3.03 x 10"^ which differs from its 
actual value by about 3%. For 7= 1.5 x 10^"*, the horseshoe 
drag is altered by the change of the vortensity gradient and 
of the density at the planet location. Again, if the difference 
between laminar and turbulent torques arose from the density 



time-evolution with turbulence, the total torque with turbu- 
lence would be w 0.85 X (AFLr + 0.55F}js) = -3.27 x 10"^ 
which is only « 8% smaller than its actual value. Note that, 
for 7= 1.5 X 10^^, the reductions of the (positive) horse- 
shoe drag and of the (negative) Lindblad torque almost com- 
pensate, which conspires to make the running-time averaged 
torque a remarkably stationary quantity. From the above 
comparison, we conclude that the torque differences at high- 
turbulence between turbulent and laminar runs can be ac- 
counted for by the time-evolution in turbulent runs of the sur- 
face density profile around the planet location. 

The time-evolution of the surface density profile with tur- 
bulence signifies that the initial smooth r^^l^ profile is not 
the equilibrium density profile for our turbulence model. The 
equilibrium profile departs significantly from a simple power- 
law profile, as suggested by the flat density transition around 
r = 0.8 at high-turbulence. However, since we observe a ten- 
dency for the time-averaged density profile to slowly steepen 
around the planet location, one may wonder whether the den- 
sity profile could tend toward a r^^l'^ profile, which would 
yield a uniform vortensity equilibrium profile. To verify this 
trend, we performed several calculations with an initial r^^/^ 
surface density profile, at high-turbulence (7 > 3 x 10^^, cor- 
responding to (ao) > 1.4 X 10"^). After about 3000 orbits, 
we find this time that the time-averaged density profile signif- 
icantly flattens out. It is possible that the equilibrium density 
profile is between the r^^l'^ and r^^/^ profiles, and we have 
not reached yet this equilibrium. In any case, it is likely that 
turbulence will tend to structure the averaged density profile, 
due to local variations in the turbulent stress. A numerical 
artifact due to our set of boundary conditions also cannot be 
ruled out. Recall that in all of our calculations, wave-damping 
zones are used along the inner and outer edges of the disk, to 
avoid reflections. This set of boundary conditions leads to 
slowly depleting the disk mass outside of the damping zones, 
in a timescale that depends on the turbulence strength. It ex- 
plains why, with increasing 7, the time-averaged density pro- 
file with turbulence takes smaller values around the planet. 
For the highest values of 7 that we considered, we found that 
forcing all fields in the damping zones toward their initial 
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Fig. 10. — Relative perturbation of the disk vortensity (with respect to its initial profile) for a laminar and a turbulent runs. Both calculations have ^ = 10 , 
hp = 5%, and = 0. Contours are shown in a narrow ring around the planet location. The left panel displays the instantaneous perturbed vortensity at 500 
orbits for the laminar run, for which a = 2 x 10^* (v = 5 x 10~^). Streamlines are overplotted as solid curves. The middle panel shows the instantaneous 
perturbed vortensity at 500 orbits for the turbulent run, for which (o!d) ^ 2 x 10^* (7 = 6 x 10"'). Instantaneous streamlines are also depicted as solid curves. 
The right panel plots the perturbed vortensity of the turbulent run, time-averaged between 400 and 600 orbits. Streamlines averaged over the same time period 
are overplotted. They are directly comparable to those of the left panel. In all panels, the vertical dashed line depicts the planet's orbital radius. 



value, as in lde Val-Borro et aTl ( 1200 6). or toward their instan- 
taneous axisymmetric value, has little impact on the density 
profile around the planet. 

Similar results were obtained with a planet mass twice as 
large {q — 10"^), as shown in the right panel of Figure [8] 
Note that these runs were performed over 2000 orbits, and 
not over 4000 orbits as in previous series with ^ = 5 x 10"^. 
At a given value of y, differences between laminar and tur- 
bulent density profiles are thus smaller with q = 10"^ than 
with ^ = 5 X 10"^, and so are the total torques. In addition, 
two difficulties arise with q ~ 10^^. First, the gap clearance 
occurs much faster. For instance, the fully saturated torque 
of the laminar run with the smallest viscosity (a « 2 x 10^^) 
increases by k, 10% between 400 and 1000 orbits. This in- 
creasing rate is about 10 times larger than for an equivalent 
laminar run with ^ = 5 x 10"^. For the smallest viscosities, 
for laminar and turbulent runs, torques were thus evaluated 
before the systematic increase due to the gap opening. Fur- 
thermore, the torque of the laminar inviscid run experiences 
large-amplitude, fast oscillations after ~ 600 orbits. These 
oscillations are triggered by the format ion of vortices flow- 
ing along the edges of the planet gap dLi et al.ll2005h . The 
amplitude of the oscillations is several times larger than the 
differential Lindblad torque, and their period is comparable 
to the planet's orbital period. We do not get any vortices and 
torque oscillations in the laminar run with the smallest viscos- 
ity (a « 2 X 

4.3. Structure of the horseshoe region with turbulence 

We have shown in ji l4.2l that. when it is time-averaged over 
a sufficiently long time period, the torque evaluated in a tur- 
bulent model is in very good agreement with the torque ob- 
tained with a similar laminar model, providing both models 
have same vortensity's diffusion coefficients (a w (ao) with 
previous notations). Under these precautions, the mean sat- 



uration levels of the horseshoe drag in turbulent and laminar 
calculations are therefore very close. This agreement indi- 
cates that, in time-average, the properties of the horseshoe re- 
gion (width, vortensity advection-diffusion) should be similar 
in both cases. To investigate these properties, we performed a 
laminar and a turbulent calculations with q = 10"^ and (7 = 0. 
The disk aspect ratio at the planet location was increased to 
hp — 5%, so that the strong vortensity perturbations induced 
by the planet wake are located outside of the horseshoe region. 
The laminar run has a kinematic viscosity v = 5 x 10^^ (in 
code units), which is equivalent to a = 2 x 10^* at the planet 
location. Similarly, the turbulent simulation was performed 
with 7 = 6 X 10"^, which corresponds to (ao) w 2 x lO""* 
from Equation dSJ. This viscosity value is low enough so that 
the surface density profiles of the turbulent and laminar runs 
are not significantly altered over the duration of these runs 
(600 orbits). 

We display in the left panel of Figure[TO]the relative pertur- 
bation of the vortensity field (with respect to its initial profile) 
for the laminar calculation, at 500 orbits. Overplotted as solid 
curves are streamlines in the planet's frame. The two par- 
ticular streamlines passing very close to the planet location 
(r — rp, (p — (pp) are the horseshoe separatrices. The vertical 
dashed line shows the location of the planet's corotation ra- 
dius rc, where D.{rc) — Clp. It also corresponds to the planet's 
orbital radius rp, as the unperturbed pressure profile is uni- 
form for these runs. Since the unperturbed vortensity profile 
decreases with radius, vortensity advection-diffusion yields 
negative vortensity perturbations along inward downstream 
streamlines ((p > (pp, r < rc), and positive vortensity perturba- 
tions along outward downstream streamlines {(p < (pp, r > r^). 
These vortensity perturbations decrease in a timescale com- 
parable to the viscous diffusion time across the horseshoe re- 
gion. As can be seen in the left panel of Figure [TOl vorten- 
sity perturbations cancel out over a time ^ Tiib/4 proceeding 
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horseshoe U-turns. 

The middle panel of Figure [TO] shows the instantaneous 
vortensity perturbation for the turbulent run, at the same time. 
For the turbulence amplitude taken here, the vortensity pertur- 
bations due to turbulence are typically larger than those ob- 
tained in the laminar run. Instantaneous streamlines are also 
depicted as solid curves, which highlight a random walk mo- 
tion rather than well-defined circulating and librating motions 
around the (expected) horseshoe region. 

However, with turbulence, we are primarily interested in 
the time-averaged fields. We thus display in the right panel of 
Figure[TO]the perturbed vortensity time-averaged between 400 
and 600 orbits. To calculate this time-average, fields outputs 
were produced every 1 /20* of orbit, which is approximately 
half the lifetime of the turbulent mode with m = 6. Solid lines 
depict streamlines averaged over the same time period, us- 
ing the same time sample. The two streamlines passing near 
the planet location show the time-averaged separatrices of the 
mean horseshoe region. All other streamlines of the left and 
right panels of Figure [TO] are calculated with the same initial 
coordinates in the r—(p plane, and are therefore directly com- 
parable. We comment that the comparison with the instan- 
taneous fields and streamlines of the laminar run is fair, as a 
steady state is reached well before 400 orbits for this run. In 
the turbulent run, fluid elements have, in time-average, librat- 
ing streamlines inside of a mean horseshoe region, and circu- 
lar streamlines outside of it. The averaged streamlines of the 
turbulent run much resemble the instantaneous streamlines of 
the laminar run. In particular, the half-width of the mean 
horseshoe region with turbulence does not significantly differ 
from that without turbulence. The time-averaged vortensity 
perturbations are also very similar to those of the laminar run. 
This qualitative agreement signifies that both simulations do 
have similar values of the vortensity's diffusion coefficient, as 
a priori expected from our relationship between y and (ao) ■ 

The above comparison leads us to the following comments. 
In time-average, turbulence tends to "diffuse" vortensity, anal- 
ogous to the effect of viscosity in laminar disks. This sim- 
ilarity could be surprising at first glance since the turbulent 
potential cannot act as a source term in the vortensity equa- 
tion. With our turbulent potential, vortensity is therefore con- 
served along instantaneous streamlines. Nonetheless, random 
motions triggered by turbulence cause turbulent diffusion, so 
that the vortensity alon g tim e-averaged streamlines is not con- 
served. As shown in § 13.41 this process can be modeled by a 
simple diffusion law, featuring the diffusion coefficient («£)). 
In the appendix, we derive the time-averaged vortensity equa- 
tion for a two-dimensional disk subject to the turbulent poten- 
tial. 

4.4. Numerical issues 

All the calculations with an embedded protoplanet were ob- 
tained with a grid resolution of 512 x 1536 and included the 
m> 6 modes truncation in the expression of the turbulent po- 
tential. We investigate here how these assumptions affect our 
results. For this purpose, we consider the turbulent model of 
§ O with ^ = 5 X 10-^ C7 = 0.5 and 7 = 4 X lO^^. Three 
additional calculations were performed, (i) one at the same 
resolution but without the m > 6 modes truncation, (ii) an- 
other one at the same resolution with cut-off of both m = 1 
and m> 6 modes, and (iii) one with the m> 6 modes trunca- 
tion but with double resolution (1024 x 3072). The running- 
time averaged torques obtained with these four calculations 
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Fig. 1 1 . — Converg ence properties of the running-time averaged torque for 
the turbulent runs of ii l4.4l The torque obtained with our standard resolution 
(512 X 1536) and with the m > 6 modes cut-off is depicted with a dotted line. 
The dot-dot-dot-dashed curve displays the torque obtained with the m < 2 
and m > 6 modes cut-off, at the same resolution. The dashed curve shows 
the torque without the m > 6 modes truncation, at the same resolution. The 
solid curve displays the torque with the m > 6 modes truncation, at higher 
resolution (1024 x 3072). 

are displayed in Figure [TT| Doubling the resolution in each 
direction has no significant impact on the torque evaluation. 
Although not shown here, we checked that this result holds 
without the m > 6 modes truncation. This convergence in res- 
olution is not surprising because our model does not include 
any feedback from the smallest length scales, where energy 
dissipates, to the largest ones, through which exchange of an- 
gular momentum primarily occurs. The reason for this is that 
turbulent modes are always regenerated independently of all 
other modes still at work in the disk. Including an energy 
equation should not alter the convergence in resolution. 

The torque r.t.a. obtained with the additional m — 1 modes 
cut-off is more negative, as expected. Discarding m — 1 
modes decreases indeed the turbulence strength, and therefore 
the equivalent alpha viscosity (ao) associated to turbulence. 
For our value of 7, decreasing (ao) weakens the horseshoe 
drag (see figure[8l), and the total torque is more negative. The 
relative difference of the total torques is « 10%. From the left 
panel of figure |8] we point out that decreasing 7 by about one 
order-of-magnitude would make the torque obtained with the 
m> 6 modes truncation almost coincide with the differential 
Lindblad torque, which is about 50 % lar ger This comparison 
underscores that, as anticipated in §[33] the absence of m = 1 
modes in the study of LSA04 cannot account for the order- 
of-magnitude difference in the expressions for the turbulent 
torque rt.a. given by LSA04 and by Equation ( fTOl ). 

Discarding the m > 6 modes cut-off renders the stationary 
rt.a. torque more positive. This is also an expected result, 
since including m > 6 modes is fou nd to increase (ur) by a 
factor of « 1.5, as shown in 13.31 Assuming that {ud) is 
increased by the same factor, and because 70c (ao)'/^ we 
expect that the r.t.a. torques obtained (i) with 7 = 4 x 10^^ 
and without truncation, and (ii) with 7 = 5 x 10^^ but with 
truncation, should be approximately the same. For the for- 
mer run, the steady torque is ^ —3.4 x 10^^ (see figure [TTli. 
whereas for the latter run, it amounts to ^ —3.35 x 10^^ 
(see left panel of Figure [8]l. This close agreement confirms 
that, for 7 = 4 X 10^^, the m > 6 modes cut-off sHghtly de- 



14 



Baruteau & Lin 



o 
a. 







w/o m > 6 cut-off 

with m > 6 cut-off 


5.5-1 0" 






5.0-1 0" 


- Y=10-''\ 


^-^7"^ Y=5x10'^ 


4.5-10" 


LY = 1.5x1 0'" V - 




4.0-10" 
3.5-10" 


: M, = 1.6M..„, 
- h = 3% 

dIogLo /dlog r = -1/2 





0.7 0.8 0.9 1.0 1.1 1.2 1.3 
Radius 



Fig. 12. — Surface density profile time-averaged over 4000 orbits for tur- 
bulent runs with and without the m > 6 modes cut-off on the expression for 
the turbulent potential. 

creases the turbulence strength. The surface density profiles 
obtained in previous cases (i) and (ii), and time-averaged over 
4000 orbits, are displayed in Figure [12] Also, we have in- 
vestigated the impact of the m > 6 modes cut-off with y = 
lO^'*. The r.t.a. torque over 4000 orbits reaches a steady 
value of ^ —1.5 X 10^^, which is approximately a factor of 
2 and 2.5 smaller than the torques with cut-off obtained for 
7 — 10^^ and 7 = 1.5 x lO^'*, respectively (see left panel of 
Figure [8]i. Clearly, for 7 — 10^^, the impact of the m > 6 
modes cut-off on the torque cannot be explained by a slight 
decrease of 7. Figure [12] shows that the corresponding den- 
sity profile (solid gray curve) significantly differs from the 
density profiles with no cut-off, and obtained for 7 — 10^^ 
and 7= 1.5 X 10^^ (dash-dotted gray curves). Interestingly, 
for this turbulence amplitude, the time-evolution of the disk 
density profile changes quite substantially without the m > 6 
modes cut-off. Around the planet location, the averaged den- 
sity profile is flatter without cut-off, which increases the pos- 
itive horseshoe drag, and it takes smaller values, which also 
makes the total torque more positive. 

5. DISCUSSION AND CONCLUSIONS 

We have revisited the interaction of a low-mass planet with 
an isothermal turbulent disk, with a special emphasis on the 
horseshoe drag desaturation by turbulence. Two-dimensional 
hydrodynamic simulations were performed, using the turbu- 
lence model of Laughlin et al. (2004). This model is based on 
applying a turbulent potential to the disk, which corresponds 
to the superposition of simultaneous wave-like modes. 

An in-depth analysis of the turbulent potential properties is 
undertaken in § [3] We show that these properties can be se- 
lected so that perturbations generated by the turbulent poten- 
tial much resemble those obtained in 3D MH D calculations. 
For in stance, the modes lifetimes taken by iLaughlin et al.l 
(120041) are reduced by a factor of ten to attain an autocorre- 
lation time of order one planet's orbital period. We quantify 
the transport of angular momentum by measuring Reynolds 
alpha parameters (ur), which we relate to the turbulence am- 
plitude 7 in Equation ((6]l. The effect of turbulence diffusion 
is also quantified by evaluating vortensity's diffusion coef- 
ficients, which we model by an equivalent alpha viscosity 
(ao) ■ A simple relation between 7 and (ao) is given by Equa- 
tion ([8]l. An important result is that vortensity diffusion is 



about four times more efficient than transport of angular mo- 
mentum: {ud) « 4(a/j). We argue in 13.41 that (ao) can 
therefore model the total alpha parameter of typical 3D MHD 
calculations (assuming they have a radial Schmidt number of 
order unity). Should our results be checked against such cal- 
culations, one would have to set up the initial conditions such 
that the average Reynolds alpha parameter (ur) be similar to 
one of our values (but note that torque differences will natu- 
rally arise from the transition from 2D to 3D). We also provide 
in Equation JTOl i an expression for the mean deviation of the 
turbulent torque distribution that gives good agreement with 
3D MHD simulations of a disk fully invaded by the magne- 
torotational instability. Our turbulence model should there- 
fore be well-suited for studying disk-planet interactions not 
only in dead zones, where angular momentum transport is 
mostly transporte d by density waves (Fleming & Stone 2003]; 
'Oishi et al."2007), but also in active layers, where turbulence 
primarily originates. Note however that our turbulence model 
has compressible modes only (turbulence is driven by a scalar 
potential), whereas MRl turbulence also contains incompress- 
ible (vortical) modes. 

By using accurate estimates of the running-time averaged 
turbulent torque, and of the fully unsaturated torque expected 
in laminar disks, we give in Equation ( i20] i an upper esti- 
mate (see § I3.6I 1 of the convergence time in our simula- 
tions. This convergence time, which is the timescale such 
that the running-time averaged stochastic torque becomes a 
small fraction of the total laminar torque, is proportional to 

— 1/2 9 

Mp ' hp. It is thus particularly sensitive to the disk scale 
height. For our conservative disk and planet parameters 
(Mp w 1 .6M©, hp = 3%,a = 0.5, t = 0), up to - 2000 planet 
orbital periods are required before the running-time averaged 
turbulent torque attains ^ 10% of the fully unsaturated torque 
expected in a similar laminar model. Increasing the planet 
mass to Mp — lOM©, and the disk aspect ratio to hp = 7%, 
yields a maximum convergence time of ^ 4500 Torb- 

We then present in §|4]our calculation results with a low- 
mass planet embedded in a turbulent disk. For comparison, 
similar laminar calculations are performed with a kinematic 
viscosity v = ac^H. The main results are the following: 

• The averaged differential Lindblad torque with turbu- 
lence takes very similar values than in laminar disk 
models, providing (i) the former is time-averaged over 
a sufficiently long time period, and (ii) turbulence does 
not significantly alter the time-averaged density profile. 

• Turbulence can unsaturate the horseshoe drag, depend- 
ing on the turbulence strength. The horseshoe drag de- 
saturation by turbulence can be modeled by vortensity 
diffusion across the time-averaged horseshoe region, 
with a diffusion coefficient D = {ao)csH. We comment 
that it is unclear how the horseshoe drag behaves when 
the largest size of the turbulence eddies becomes com- 
parable to, or larger than the width of the horseshoe re- 
gion. Vortensity should then enter the mean horseshoe 
region in an advective way, rather than in a diffusive 
way. This situation deserves attention, particularly for 
low-mass planets. It is possible that the horseshoe drag 
value is determined by the advection timescale across 
the mean horseshoe region, at the averaged turbulent 
velocity, in comparison with the libration and U-turn 
timescales. 

• For similar vortensity's diffusion coefficients (a w 
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(«£))), time-averaged total torques with turbu lence are 
compared with steady laminar torques in § 14.21 (see 
Figure |8]i. At low-turbulence, turbulent and lami- 
nar torques are in very good agreement. At high- 
turbulence, differences arise, which can be fully ac- 
counted for by the time-evolution of the averaged den- 
sity profile with turbulence. These results have two 
implications. On the one hand, the torque felt by a 
low-mass planet can be decomposed into a laminar 
torque, and a stochastic torque due to turbulence. On 
the other hand, the averaged value of the horseshoe 
drag in turbulent and laminar disk models are very sim- 
ilar. The same comment applies to the structure of the 
mean horseshoe region, as shown in § 14.31 It indi- 
cates that the criterion for the horseshoe drag desatura- 
tion, Tu_tum < Tvisc < ^ib, Still holds in turbulent disks, 
in time-average, at least for the turbulence amplitudes 
considered in our study (but see the above comment on 
the advective regime). The quantity Tvisc corresponds 
to Tvisc ~ /D, where D — {aD)csH is the vortensity's 
averaged diffusion coefficient, and is the half-width 
of the mean horseshoe region. 

With turbulence properties as close as possible to those of 
3D MHD simulations, we find that the horseshoe drag exerted 



by isothermal disks on low-mass planets can remain unsatu- 
rated on the long term, depending on the turbulence strength. 
These results require confirmation by 3D MHD long-term 
simulations of planet-disk interactions, with disks either fully 
magnetized or harboring a dead zone. It would also be of 
relevant interest to investigate how accretion onto a low-mass 
planet can be affected by turbulent motions in the planet vicin- 
ity. In forthcoming works, we will extend our study to radia- 
tive disks, and we will investigate the trapping of a protoplanet 
at a density transition in presence of turbulence. 
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APPENDIX 

TIME- AVERAGED VORTENSITY EQUATION 

We derive in this section the time-averaged vortensity equation for a two-dimensional disk subject to the turbulent potential 
described in §[3] In a frame rotating uniformly at angular velocity Q., the continuity and momentum equations are 



dt 



EV-v==0, 



— + (vV)v- 



-2nkxv=-V(/i + <I>), 



(Al) 



(A2) 



where v is the two-dimensional velocity field. In Equation (IA2b . h = J dp /T, is the fluid enthalpy, <I> is the total potential felt 
by the disk (including the time-dependent turbulent potential), and k is the unit vector in the vertical direction. Using the vector 
identity (v • V)v = (1 /2) Vv^ — v x w, where w = V x v is the vorticity, and taking the curl of Equation (IA2l i. we find 



(9(w + 2ak) 
di 



-Vx [(w + 2nk) xv] = 0. 



(A3) 



(A4) 



Multiplying Equation ( |A3t by E \ and subtracting the product of (w + 2i2k) /E^ and Equation ( lAll l. we are left with 

| + (vV)^^0, 

where ^ = (w + 2i2k) /E is the vorten sity. In the above derivation , we used the vector identity V x (A x B) = A • ( V • B) — (A • 
V)B — B(V-A) + (BV)A. Equation (IA4l i reduces to equation 2 in Koryc ansky & P apaloizou ( 1 996) in the case of a steady flow. 
We use the so-called Reynolds decomposition, in which an instantaneous value is written as the sum of a mean value (denoted 
with a zero subscript) plus a fluctuation: v = Vo + 5v and ^ = + 5^. Mean values are taken over a timescale T that is large 
compared to the correlation timescale, but short compared to that of the flow evolution. We thus have v : 
d\ = 0. Similarly, £, — £,o and 5£, = 0. After time-averaging Equation ( IA4I ). we find 

§ + (..v)T 



T ' Jq \dt = vo and 



-(5v-V)5^ =0, 



(A5) 



where the time interval between two successive values of t is greater or equal than T, and where the term (5v- V)5i^, which 
features the correlation between 5y and 5^ , is responsible for the non-conservation of vortensity along time-averaged streamlines. 
It is analogous to the Reynolds stress in the averaged momentum equation. 
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